Dynamic hyperpolarizability of the one-dimensional hydrogen atom with a 5-function interaction 
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The dynamic hyperpolarizability of a particle bound by the one-dimensional 5-function potential is obtained 
in closed form. On the first step, we analyze the singular structure of the non-linear response function as given 
by the sum-over-state expression. We express its poles and residues in terms of the wave-number k. On the 
second step, we calculated the frequency dependence of the response function by integration over k. Our method 
provides a unique opportunity to check the convergence of numerical methods, and is in a perfect agreement with 
the static and high frequency limits obtained by different theories. The former is obtained using the approach 
of Swenson and Danforth (J. Chem. Phys. 57, 1734 (1972)). The asymptotic decay is studied using the method 
of Scandolo and Bassani (Phys. Rev. B 51, 6925 (1995)). Its extension to the case of quadrupole polarizability 
reveals a universal (not dependent on the choice of the system) asymptotic behavior of the hyperpolarizability. 

PACS numbers: 42.65.An 



I. INTRODUCTION 

There is only a small number of exactly solvable realistic 
quantum systems. The hydrogen atom is one of the most ex- 
perimentally and theoretically studied. Remarkably, expres- 
sions for the Green function [1], transition matrix elements 
between the ground and excited states [2], and a number of 
other properties of this system exist in closed form. How- 
ever, the complexity of the matrix elements as well as of 
the sum-over-states (SOS) expressions for the nonlinear op- 
tical response functions (3J hinders an analytic computation 
of its dynamic hyperpolarizabilities. Already expressions for 
the dynamic polarizability of hydrogen atom are very compli- 
cated Hi. 

To describe the multiphoton ionization of atomic hydrogen 
it is sufficient to consider summations over all intermediate 
bound states [5fl. In contrast, the computation of the nonlin- 
ear optical response is much more involved since summations 
over the continuum states should also be performed |0] . This 
considerably complicates derivations and calls for the devel- 
opment of new methods. Thus, Mizuno [7] used the sturmian 
coulomb Green function (8J[ in order to compute the third har- 
monic generation coefficient of the hydrogen atom. These re- 
sults were extended by Shelton [9] to other third order pro- 
cesses. Although written in analytical form, they still con- 
tain infinite sums, and thus can only be analyzed numerically. 
Therefore, it is desirable to have a simpler model system that 
would allow for a closed-form solution for non-linear optical 
responses. 

A particle bound by a one-dimensional 5-function potential 
bares a close resemblance to the three-dimensional hydrogen 
atom with Coulomb potential. Despite its simplicity, the sys- 
tem is attractive for the theoretical analysis because like the 
hydrogen atom it contains both discrete and continuum eigen- 
states. Also the Green function of these two systems shows 
remarkable analogies [10]. The (5-potential allows for a great 
simplification of the analytic work since one only needs to 
carry out the integration of relatively simple functions. 

The frequency-dependent electric polarizability of this sys- 
tem has previously been computed by explicit solution of 



the time-dependent Schrodinger equation with uniform time- 
dependent electric field and expansion of the charge density in 
the perturbed ground state up to and including terms linear in 
applied electric fields [11]. Recently, the result has been con- 
firmed using a perturbation technique 11211 based on the work 
by Nozieres 11311 . Here, we present results for the lowest order 
non-zero hyperpolarizability due to the combination of two 
electric-dipole and one electric -quadrupole transitions. This 
response function describes the second-harmonic generation 
(SHG) process. SHG due to only dipole transitions is forbid- 
den because the system possesses an inversion symmetry. The 
next order non-linear response (the third harmonic generation) 
due to the combination of four dipole transition should have a 
comparable magnitude and will be analyzed elsewhere. 

As a starting point, in Sec. II we study singular properties of 
the SHG response function as given by equations of Orr and 
Ward H. We assume that all conditions of validity of these 
formulas are fulfilled. The reader is referred to the book of 
Shen[14] for the discussion of relevance of these expressions 
for a certain experimental situation. Using expressions for 
the poles and residues of the response function obtained here 
we perform the actual computation for the one-dimensional 
hydrogen in Sec[TTT] Finally, we qualitatively analyze obtained 
results in the low and high frequency limits and compare them 
with prediction of other theories (Sec. lIVb . 

In this respect we refer to the work of Fernandez and Cas- 
tro 11511 which gives an introduction to the method that avoids 
an integration over the continuum spectrum. Based on the 
original idea of Swenson and Danforth 11611 . and extending the 
derivation of Austin [17], they obtained a recursive equation 
that relates the expectation values of x N operator at different 
level of the perturbation expansion. They used the method to 
find a perturbative correction to the ground state energy and 
to describe the appearance of resonant states in the presence 
of static electric field. As will be shown below, the same ap- 
proach can be used to obtain an arbitrary order static hyper- 
polarizability 11811 . We performed corresponding calculations 
for our model system and found an agreement with the static 
limit of our SOS-derived expression. 

The high frequency limit of the second-harmonic genera- 



2 



tion in dipole approximation was obtained by Scandolo and 
Bassani 11911 . Later, this result was generalized to higher har- 
monic generation processes using the general quantum theory 
of Kubo optical response function [ 20] . It yields l/ui 2n+2 
decay of n-th order response function. We show that non- 
local response (the quadrupole polarization) dictates a differ- 
ent asymptotic behavior, 1/ui 4 for the SHG process, in agree- 
ment with our SOS results. 

In AppendixlAl we outline calculations of the transition ma- 
trix elements for this system, and in Appendix iBl we derive 
the frequency dependence of the SHG response function from 
the corresponding residues. An extensive use of the complex 
analysis and generalized functions is made. All calculations 
are done analytically. The final results are expressed in terms 
of elementary functions. 



II. SINGULARITIES OF THE SHG RESPONSE 
FUNCTION 

The poles and residues of the second-order response func- 
tions can be derived from the microscopic expressions for the 



nonlinear polarizations (as, for example, given by Eq. (43) of 
Ref. [H) assuming certain form of the light-matter interaction. 
It is important to notice that corresponding perturbation op- 
erators V as well the polarization operator P must be used 
in the renormalized form with vanishing expectation values in 
the ground state 
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where V and P are the bare operators. Other forms of pre- 
senting the non-linear response functions are also known. For 
example, lifting the requirement on the diagonal elements 
to be zero leads to the appearance of additional terms l2lll . 
They can be traced back to the secular divergence problem. 
Clearly, both forms lead to identical results at the end. To 
make our manuscript self-contained, we reproduce here the 
equation of Orr and Ward for a general second-order process 
(u> a = u>i + L02) using slightly simplified notations: 
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The simplification concerns the use of only renormalized op- 
erators, while in the original formulation the authors had both 
renormalized and bare operators. Such transformation is al- 
ways possible since (V) ng = (V) ng (the same also holds for 
P). As a consequence, the expression acquires a more sym- 
metrical form. Here, Jx,2 denotes the average of all terms gen- 
erated by permuting w\ and u>2- K is a numerical factor that 
depends on the permutational symmetry of incident photons. 
For the SHG process it is equal to h. The summations are run- 
ning over all excited states excluding the ground state. 9 mg 
denotes the energy difference between the excited state m and 
the ground state (labeled by g). In order to avoid divergencies 
at the resonances, a phenomenological damping (iT) of ex- 
cited states is introduced, Q mg = E m — E g + iT. This shifts 
the poles of the response functions away from the real axis on 
the complex plane. 

If the damping constant iT is small, it has no relevance for 
the derivation of residues. That is why we provisionally as- 
sume that all poles lie on the real axis (f2* g = O nff ). Fur- 
thermore, a simple analysis of Eq. OJ shows that the poles 
of nonlinear polarization are symmetrically situated around 
u) = 0. Additional symmetry also exists for their residues. 



The residue of a pole at negative energy has opposite sign to 
the residue of its positive energy counterpart. Thus, for our 
analysis it is sufficient to consider singularities at positive en- 
ergies only. 

A simple analysis of the SHG case (u>i = u>2 = uj, w a — 
2w) shows that each excited state m gives origin to 4 poles 
on the complex plane uj = ±f2 mfl and u> = ±Q mg /2. The 
SHG response function is defined as a second derivative with 
respect to the external electric field 
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Its residues for positive frequencies will be denoted as fol- 
lows: 



4TV) = Res w=1/2fimsX gr ) . (5b) 



There will be two terms contributing to the residue of the po- 
larization at u> = n ng : 
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Here, we substituted the frequencies and the symmetry factor 
K and expanded the symmetrization operator. A bar over the 
matrix elements denotes their symmetrization, which signifies 
the equivalence of two incident photons and follows from the 
application of 1\ 2 to the whole expression. In the same way, 
we write an expression for the residue of the polarization at 
u> = Vt mg /2 where only one term contributes: 
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Equations © and (0) are obtained from the general the- 
ory, and, therefore, are valid for systems of arbitrary di- 
mensionality and for different light-matter interaction mech- 
anisms. They are not only useful in present theoretical anal- 
ysis, but also can bring a substantial computational savings 
when used in ab initio calculations for realistic systems. This 
directly follows from the estimates on the number of floating- 
point operations needed to directly evaluate Eq. [3] (N sos — 
0{N U1 ■ N 2 )) in comparison with a two-step procedure, where 
the residues (Eqs. \6\1\ are computed with numerical cost of 
Nl os = 0(N 2 ) only, and the frequency dependence is ob- 
tained on the second step using N** a = 0(N UJ ■ N) opera- 
tions. Here we assumed that the system has N excited states 
and N u frequency points are required. 

On the last step, we specify Eqs. ( 1617b for the model system. 
Since the SHG is strictly forbidden in centrosymmetric sys- 
tems within the electric dipole approximation (Xeee = 0), 
we will be considering a one-component response function 
Xqie w i tn the dipole perturbation operator V — —exF, and 
the quadrupole polarization P — ex 2 . Furthermore, we make 
use of atomic units, the conversion is done by setting the elec- 
tron charge e = — 1, and m = K = 1. Finally, the differentia- 
tions with respect to the external electric field F are performed 
according to the definition of the response function (Eq. |4j». 

We will denote two contributions to A^^\m) as A I (k) 

and A n (k), and B^'j?\m) will be denoted as B{k). It is 
natural to name them as spectral functions of the second-order 
nonlinear response in analogy with the spectral functions in 
the many-body perturbation theory. 
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Instead of a discrete state number m, these functions now de- 
pend on the wave-number k > that characterizes excited 
states of the system as will be shown below. In order to be 
consistent with naming of states of the one-dimensional hy- 
drogen that will be introduced below we changed the notation 
for the ground state to 0. We also introduced the notation 
Aa; 2 = x 2 — (x 2 )qo and took into account that (a;)oo = for 
our model system. 

III. DYNAMIC HYPERPOLARIZABILITY OF THE 
ONE-DIMENSIONAL HYDROGEN 



The system discussed here is described by the one- 
dimensional time-independent Schrodinger equation that can 
be written in the dimensionless form: 



H o = -4>"/2 - 5{x) 



Ecj). 



(9) 



The solution of the eigenvalue problem (Eq. [9j» yields a sin- 
gle bound state with the energy and wave-function: 



So = -1/2, 0o(aO=exp(-|aT|), 



(10) 



and a continuum of unbound states |[ITl,l22tl . Since the Hamil- 
tonian of the system is invariant with respect of space inver- 
sion, the wave-functions can also be constructed to have a well 
defined parity: 
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They are degenerate and have the energy E(k) = k 2 /2 in re- 
semblance with the unperturbed states of a free particle. How- 
ever, in contrast to the free-particle case the wave-number k 
can only attain positive real values. One can readily demon- 
strate the completeness and the normalization of the above set 
of eigenfunctions. 

As a starting point for A 1 (k), A 11 (k) and B(k) computa- 
tion, we need the bound-free and free-free transition matrix 
elements of the electric dipole: 
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Vl + fc' 2 

For the non-local SHG, we also need the quadrupole transition matrix elements between the ground and even unbound states: 
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and between two free odd states: 
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The derivation of the matrix elements is done in AppendixlAl 

The spectral functions are computed according to Eqs. ([8j, replacing the summations with integrations over the wave-number 

(fc'): 
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where we introduced notation uo(k) = E(k) — E = (fc 2 + l)/2. Integrals Ai(fc) and A 2 (fc) are computed in AppendixlAl 
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Here, we use twice the integration by parts and the fact that the integrand vanishes at the ends of interval. Thus, the final 
expression for A(k) — A 1 (fc) + A II (k) results from the addition of Eq. (fT~5T > and Eq. ( TTol i: 
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Computation of B(k) is slightly more involved since the resulting function appears to be non-analytic: 
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FIG. 1 : SHG spectral functions of the one-dimensional hydrogen with a 5-function interaction. All quantities are given in atomic units. Three 
panels shows three possible excitation pathways leading to the quadrupole polarization of the system. On the insets, thick lines denote resonant 
transitions. Dashed lines correspond to 2ui transitions. Eq = —1/2 denotes the ground state of the system. Unbound states have parabolic 
dispersion. 



Before the evaluation of Xqef it is instructive to analyze the properties of the spectral functions A I (k), A 11 (k), and B(k). 
They describe different excitation pathways (Fig.Q]) in the system. A peak slightly below k = \ related to the existence of the 
excitation threshold is common for them. It is interesting to observe that A 1 (k) and A 11 (k) decay as k~ 8 for large wave-number 
k. There is, however, a cancellation of two terms in their sum leading to the same asymptotic behavior as B(k) (fc~ 10 ). 

We compute the SHG response using the following representation: 
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(19) 



where T > is infmitesimally small broadening of the states that ensures the causality of the response function. We transform 
the integration to the whole real axis using the fact that A(k) and B(k) are even functions. In this form the residue theorem can 
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FIG. 2: SHG response functions of the one-dimensional hydrogen with a 5-function interaction. The axes are labeled in atomic units. Dotted 
curves represent parabolic dispersion of the response functions in the vicinity of uj — 0, as given by Eqs. l !25l26t . 



easily be applied. It is convenient to introduce following auxiliary functions: 
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and to split the integration of B(k) into two parts: 
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Thus, the response function can be written as Xqee = A(lu) + B{io), with A{uj) ~ a{oj) + a*(—uj), and B(lu) = b(cu) + 
b*(—uj). Although one can apply Sokhotsky's formula in order to quickly get its imaginary part, we will explicitly perform 
the integrations in order to obtain the real part. The calculations are quite cumbersome involving numerous applications of the 
residue theorem. Since A4(fc) is not an analytic function, we were forced to use its integral representation and to exchange the 
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order of integrations. Details of these calculations are presented in AppendixlBl Thus, final expressions are given by: 

t lit \ 1 / (1 - 12MW2M - 1 U > h, , 
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Rei(w) = — -= + T7^zS VI + 2^(1 + 12w) - Vl - 2u>(l - 12w) bl < 1, (23b) 
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IV. LOW AND HIGH FREQUENCY LIMITS 

Qualitatively, the properties of obtained response function 
(Fig. I3 can be understood from general principles. As ex- 
pected, it contains both real and imaginary parts related by the 
Kramers-Kronig relations, which are well known to hold also 
for the second harmonic response [19]. The real part of the 
function is associated with the conversion between the funda- 
mental and second-harmonic frequencies, while its imaginary 
part is associated with the removal of energy from the electro- 
magnetic field 112 311 . In our case, the Kramers-Kronig relations 
directly follow from Eq. dl9V No optical transitions are possi- 
ble for the energy of photon below the gap between the ground 
and unbound excited states. Therefore, no energy absorption 
is possible in this regime. Im^4(o;) vanishes for to < 1/2 and 
InxB(u;) vanishes for u> < 1/4. ReA is sharply peaked at 
the edge of the continuum: ReA(^) = —9.33, its frequency 
derivative has an infinite discontinuity Re^4'(o; = |— ) = 
— co, Re^4'(o; = |+) = 213.37. Similar behavior is found 
forRcS atw = \: RcB(\) = 52.14, Rc£?' (a; = \-) = oo, 
ReB'(u) = j+) = —2365.37. Similarly to the imaginary 
part, the real part of the response function remains finite and 
continuous for all values of u> and has a well defined static 
limit: 
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was already found in many realistic systems and proved ana- 
lytically 0S. 

Let us now obtain the static limit closely following the 
derivation of lfl5ll . We consider the one-dimensional station- 
ary Schrodinger equation: 

Hip = eip, H = -D 2 + V(x), D = d/dx, (28) 

where the boundary conditions are supposed to be 

ip(x — > ±oo) = 0, 

and V(x) = —S(x) + Ax. Using the method of Swenson and 
Danforth ifliitl one can obtain recursive equations 
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for the expectation values of of a;", i.e. JfW = {x N ) at 
different orders of the Rayleigh-Schrodinger perturbation the- 
ory: 
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lW=^xf»y. (31) 
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From X<, 0) = 1 and X ( 1} = one obtains with a help of 



Eq. d29jl 
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720, 
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At small lo, one observes a quadratic dependence of the SHG 
response on the frequency of incident photon. This behavior 



Starting from these values other coefficients can be recursively 
computed using Eq. ( f30b . Simple symmetry consideration 
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show that X 
nonzero values are 
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if N + p is an odd number. First few 
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X^> = 744, Xf 5 = 36960, 
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(3) 



-184080. 



The knowledge of X p allows to compute the static response 
functions. Let us make a rescaling q = x/2, cf>(q) = ip(x), 
F = — 4A, E = 2e, and Eq = 2e = — h m order to bring 
Eq. (f28]) to the form: 

Ho^q) = -0"/2 - 5(q)<f>(q) - Fqftq) = E<f>(q). (32) 



As a consequence an expansion 



must be compared with Eq. ( f3Tb yielding: 



(33) 



Q P N) = (-i) p 

Hence, the first few values are: 

a(0) 



X 



<N) 



2 n 4p ' 



(34) 



2 X 



(0,0) _ ^(2) _ 
qee ~ ^2 ~ 



1 - 

r (2 



5 
4' 



2 2 4 2 



93 



We see that Q\ correctly reproduces the static polarizability 
of our model system. The frequency dependent polarizabil- 
ity was obtained in ifTll [l2ll . Taking a limit ui — > yields 
a(0) = |. There is an additional factor of 2 that must be mul- 
tiplied with Xqe? in order to equate it to Q% . Appearance of 
this factor is due to different representations of electric fields 
for static and dynamic case used in [3]. The resulting differ- 
ence in terminology is discussed in details in l26ll . alterna- 
tively one can see it from the difference in factor K(u>i, 0J2) 
adopted by Orr and Ward for static and dynamic cases (for 
SHG K(uj,uj) = 1/2, while in the static case if (0,0) = 1, 
see Tab. 1 of M). 

Following the prescription of [ 19] the uj — > 00 behavior will 

in terms of its Fourier 



be obtained by i) defining the xk 
transform: 



Aqee 



J dr+ J dT-G ( - 2) {t 1 ,t 2 y uT+ , (35) 



where r + = t\ + t 2 , T = (t± — t 2 ) /2 and 

G^(t u t 2 ) oc -/(ti,* 2 )ff(ti,i2), (36) 
f{t u t 2 ) = 6{h)e{t 2 - h) + 0(t 2 )0(ii - t 2 ) (37) 

ii) integrating Eq. (l33T > by parts 



v (">>"0 
Aqee 



1^ f-i^m+l (38 > 



(-iu) 



iii) seeking the lowest nonvanishing order of the time deriva- 
tive gf+n g{ti, t 2 ) of the correlation function 

g(ti,t 2 ) = {[x(-t 2 ), [xi-h)^ 2 }]^. (39) 

Since we consider here the one-dimensional case, spatial in- 
dices are omitted. After some algebra one can show that the 
first nonzero term is: 



Ot 



d 2 , n d 2 g 



dhdt 2 



(40) 



Thus, the first nonvanishing derivative in Eq. (138t is of the 
third order 



d 3 

a7+3 



G (2 \h,t 2 ) oc 4[«(*i)0(t a - h) + 5(t 2 )0(ti - t 2 )\ 



By inserting this expression in Eq. d38l ). performing the inte- 
gration and taking the limit we obtain 



(cj,o;) 
Xqee ^ 



1 



for lj 



00. The same behavior is seen indeed from 



Eqs. (123124b . One can also note a delicate cancellation of oj~ 2 
terms in the sum A(oj) + B(ui). 

Another important conclusion that immediately follows 
from our asymptotic analysis is the prefactor of 1/w 4 de- 
cay. In the SHG case due to dipole transitions the response 



function behaves as 



d-'v 



)o~^, where the av- 



i,jk \dxidxjdxk 1 

erage is performed in the ground state of the system [H 
The potential V describes the electron-ion interaction. For 
the third harmonic generation one obtains in the same way 



Xi,jhl 



a 4 v 



}ott» 12711 ■ It is clear that in both 



^ dxidxjdx^dxi ' 

cases the potential and its partial derivatives are system spe- 
cific. In contrast, the quadrupole SHG response function does 
not contain derivatives of the potential, it only depends on the 
system density in the high frequency limit. 



V. CONCLUSIONS 

In this work, we analyzed the singular structure of the 
second-order response function based on the sum-over-states 
expression. As an illustration of the method the non-local sec- 
ond harmonic generation is obtained for the one-dimensional 
hydrogen atom with a <5-function interaction. It provides a 
new paradigm where the sum over continuum states can be 
analytically evaluated. Finally we analyze the static and high- 
frequency limits, which we are also able to compute using 
different methods. 

Our results lead to the fast algorithm for the calculation of 
non-linear response using the sum-over-state approach. The 
analytic expression for the SHG response can be used as a 
valuable test of the numerical convergence of SOS methods. 
The analysis of the high-frequency behavior shows a universal 
oj~ 4 behavior, which can easily be tested experimentally. 
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APPENDIX A: MATRIX ELEMENTS 



One of the simplest matrix elements is the bound-free transition dipole moment. The integrand belongs to the space of 
square-integrable functions. The integration can be done using elementary methods: 



<0o|a#-(fc)> 



-= [ dxe ^xsin(kx) = — 7= -77- f dxe x coa(kx) 

V* J-oo V* dk Jo 



^{k 2 + \f 



(Al) 



Similarly, one obtains the matrix element for the quadrupole transition between the ground and even excited state: 



(MA<t>- 



1 1 f 00 / i 

= —= - = / dxe~^x 2 cos(fcx) - - sin(/c|a;| 

2 1 ( d2 I™ a ^ 1 d2 f°° a -* ■ n 

= = ■ = -—7 / dxe cos(kx) — — —-tt / dxe smikx 

V* VTTk 2 \dk 2 J kdk 2 J 



v^vTTfc^(i + fc 2 ) 2 ' 



(A2) 



All other matrix elements mentioned in this manuscript cannot be represented in terms of only elementary functions. This is 
due to the fact that free particle states are not square integrable, but rather normalized on the ((-function. Thus, the quadrupole 
transition moment between the odd states can be computed as follows: 



oo 2 °° 

(<j>-(k')\x 2 \<j>-{k)) = - f dxsm(k'x)x 2 sm(kx) = ---^^ I dx (e ik ' x - e~ ik ' x ) (e ikx 
7T J 7T (2%Y dk z J V / 



-ikx\ 



dk 2 



[6(k + k 1 ) - 5(k' - k)] = -8 {2 \k' - k), 



(A3) 



where the last transition is valid because of k, k' > in our case. The calculation of the dipole transition moments between the 
free states of different parity can be done as follows: 



(<Mfc)M<Mfc')> 



i k 

__k 



dx cos(kx)xsm(k / x) — — j dx sm{k\x\)xsm{k' x) 



S^\k'-k) 



k' 



7T {k 12 - k 2 ) 2 



using following integrals 

oo oo 

- / dx cos(kx)xsm(k'x) = -— — / dx ( e lk ' x + e ~ lk ' x ) ( 

IT J 7T 2 2 dk' J \ 



— [S(k + k') + 5(k' - k)} = Vk 1 - k), for k, k' > 
dk' 



oo oo 

— / dx sin(fc|a;|)a;siri(A:'a;) = —7-^77 / dx sgn(x) s'm(kx) cos(k / x) 
k J k dk' J 



(A4) 



1 1 d 

k Ai dk' 

l_d_ 

k~dk' 



dxsg 



n(x) (e lk ' x 



+ e 



ik x \ ( Akx 



e — e 



-00 

1 



k' + k k' - k 



-ikx\ 



4k' 



dk' k' 2 - k 2 (k' 2 -k 2 f 
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v Im q 



Im q 




A g = -i^i(l-* 2 ) 



FIG. 3: Integration of A4(w). Poles that contribute to the integral are shown in black. 

9 = ±A 



= ±i and q = ±fc are second order poles. 
i(fc 2 — 1) are simple poles that lay on the imaginary axis for k < 1 and on the real axis for k > 1. One can verify that poles of each 
pair have residues of opposite sign; therefore, their contributions to the integral cancel if the singularities are situated on the real axis. 



where we used a property dx sgu(x)e lkx — 2i/k. The spectral functions A 1 (k), A 11 (k), and B(k) contain contribution 
from four integrals (Xi (fc), i = 1 . . . 4). Integrals Ai (k) and A3 (k) are evaluated by the integration by parts and A2 (fc), A4 (fc) are 
evaluated by the contour integration along the real axis closed by an infinite semi-circle in the upper half of the complex plane 
(Fig.0, as follows: 



Ai(fc) = 



A 2 (fc) 



A 4 (fc) 



1 



q 2 dd{q - k) 



2k 2 -q 2 + l (1 + q 2 ) 3 dq 
u 

2fc(fc 2 - 1) 



d 1 = - I % - k )-j^ 



2k 2 - q 2 + 1 (1 + g 2 ) 3 



dq 



(k 2 + l) 5 ' 

X) 

1 



2k 2 -q 2 + l (1 + q 2 ) 3 (k 2 -q 2 ) 



2 dq 



1 



2 q 



1 



2k 2 - q 2 + 1 (1 + q 2 f (k 2 -q 2 ) 2 

00 

q dS(q — k) 



2k 2 - q 2 + 1 (1 + q 2 ) 3 {k 2 - q 2 ) 2 
7T fc 4 + 12k 2 - 23 



dq 



32 (fc 2 + l) 5 ' 



2<? 2 - fc 2 + 1 (q 2 + l) 2 dq 



Ik 2 -I 



-dq 



S(q-k) 



dq 



1 



2q 2 -k 2 + l{q 2 + l) 2 



dq 



(k 2 + l) 4 



1 







2q 2 -fc 2 + l(q 2 + l) 2 (q 2 -fc 2 ) 



2 d 1= 



2q 2 -k 2 + l(q 2 + l) 2 (q 2 -k 2 ) 2 



dq 



-2m 
2 



Res(fc = i) + Res ( k = i 
Res(fc = i), 



1-fe 2 



fc < 1 
fc > 1 
I — 



4?r 



7T fc 2 -7 



4 (1 + fc 2 ) 4 (1 + fc 2 ) 4 



1-fe 2 



(A5) 



(A6) 



(A7) 



, fc < 1 
fc > 1 



■(A8) 



APPENDIX B: HILBERT TRANSFORMS per complex half -plane encircled by the contour of integra- 

tion depends on the value of frequency uj (Fig. |4). When 
u> < h, i.e. the energy of incident photon is insufficient to 
Below, we evalu ate integrals q (u), bi(u) and b 2 {uj). They generate optical transition from the ground state to the con- 
are defined by Eqs. d20|22al22bt . respectively. These are quite tinuum of unbound states. The imaginary part of the function 
cumbersome calculations, which involve application of the a (w) vanishes. While the real part has a contribution from two 
residue theorem. For a(u>), the number of poles in the up- 
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Im k 



Im k 



k = 




k = s/%> - 1 Re k 
k = —i 




FIG. 4: Integration of a(u>). k = ±i are poles of the seventh order, k = ±^J2u — 1 are simple poles. It is important to notice that due to the 
presence of small imaginary part (cu — > lu + iF, F > 0) they are shifted away from the real axis even in the case of lu > i. Thus, only one 
pole of this pair is encircled by the contour for any value of lu. However, when lu > i, its residue is real, thus contributing to the imaginary 
part of the integral. In all other cases, the poles have imaginary residues. This fact is reflected by different notations for singularities with real 
(crossed circle) and imaginary (black circle) residues. 



poles. For the energy of an incident photon above the ioniza- one pole contributes to the imaginary part of the function, 
tion threshold, i.e. to > |, one pole contributes to the real and 

I 



Rea(w) = iri 



Res k=i 
Res fc=J + Res k=lVT ^ 



LU > 
LU < 



r 



A(k) 



lu — LUo(k) + iT 



U7lu 6 + IOGlu 5 + 86lu 4 + 88w 3 + 184w 2 - 208w +16 1 



Ima(w) = ni 



Res t _ 



k=y/2UJ=l w > 2 







LU < 



256w 7 

_L 

2 
1 







LU > 



V. 



A(k) 



16lu 7 \ VI - 2aJ(-12w +1) lu < 2 
1 \ y/2u- l(-12w + l) lu>L 



lu — LUa(k) + iT 16tt> 7 







LU < 



Calculations for b\ (lu) can be done along the same line with a small distinction that now lu = | is a point, in which the 
function changes its behavior. This is because b\(uS) describes resonant 2lu transitions. The result works out to be: 



Rebi(w) 
Im bi(uj) 



28lu 6 + Ulu 5 + 8lu a + 6lu 3 + IOlu 2 - lieu + 2 1 



64w 7 







64w 7 1 VI - 4w lu < i; 



64w 7 



V4w - l(7w - 2) lu > \ 



LU < 



r 



The computation of b 2 (lu) brings us to the following double integral: 



32 

&2M = -T / 



dk- 



(q 2 + l) 2 7 (1 + k 2 ) 3 (w - Wo(fc)/2 + tT/2)(2q 2 - k 2 + l)( g 2 - /c 2 ) 2 ' 

— 00 —00 

After simple but lengthy consideration of all cases (Fig. [5} we obtain: 



ts(w) 

Re6 2 M = s(lu) +u 2 (lu) + { 

ui(u) + u 3 (lj) w < \; 



5 > w > i 



+ 





LU > 



lmb 2 (u) = { h(Lu) + t 2 (Lu) ±>w>±, 
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FIG. 5: Scheme that shows how the double integration of 62 (w) is performed. In both cases, the contour of integration encircles the upper 
complex half-plane. On the first step, we find the residues of poles contributing to the integral over k. Furthermore, we consider all possible 
ramifications that arise for different values of to when integrating over q. Some branches of the integration tree lead to the same resulting 
function. In order to disentangle contributions to real (imaginary) part of the integral, we denote residues that have real (imaginary) value by a 
single (double) frames. 



where the auxiliary functions (s(co), ti(w), and i = 1, . . . , 3) are given by: 

197w 2 + 50lu + 10 



s(uj 
ti{u> 
t 2 (w 
t 3 (uj 

Ui(u> 

u 2 {uj 
u 3 (uj 



Rcs q=l (Res k =i[f(k,q,uj)]) = ;r . ( . 
Res 9=J (Rfis fc= v3s=T[/(fc,g,w)]) = — -, 



Res 9=v ^T (Res fe=v q^[/(fc,g,w)]) = ^ , 



Res g=i (Res fc=l/I 34^[/(fc, 9, w)]) 



Vi - - 2) 



Res 9=lv ^^j (Res fe=lVT ^[/(fc, g , w)]) = Res 9=V 4^ (Res fe=v q^ T [/(fc, g, w)]) 



7cj - 2 



= Res 9=iv ^^j (Res fe=lVT 34^[/(fc, 9, w)]) 



VI - 4^V1 - 2w 
1&7 • 



Here we denote the function under the integrals as f(k, q, u). 
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